Computing absolute free energies of disordered structures by molecular simulation 
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We present a Monte Carlo simulation technique by which the free energy of disordered systems can 
be computed directly. It is based on thermodynamic integration. The central idea is to construct 
an analytically solvable reference system from a configuration which is representative for the state 
of interest. The method can be applied to lattice models (e.g., the Ising model) as well as off-lattice 
molecular models. We focus mainly on the more challenging off-lattice case. We propose a Monte 
Carlo algorithm, by which the thermodynamic integration path can be sampled efficiently. At the 
examples of the hard sphere liquid and a hard disk solid with a defect we discuss several properties 
of the approach. 
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The fundamental equation S = f{U,V,{Na}), which 
connects the entropy S with the internal energy U, the 
volume V, and the numbers of particles of type a, 
contains all information about a system that is accessible 
within classical thermodynamics. Other thermodynamic 
potentials such as e.g. the free energy are related to the 
fundamental equation by Legendre transform, and hence 
they equally contain this informationi.. Therefore there 
is large interest in computing free energies in many ar- 
eas of science, i.e., statistical physics, materials science, 
theoretical chemistry, and biology-. 

There are only very few, special cases in which the free 
energy of a system can be computed directly: Either the 
accessible phase space volume can be enumerated com- 
pletely (as e.g. for a lattice gas model on a small lat- 
tice), or the problem can be solved analytically in the 
first place (as e.g. for the ideal gas). In all other cases 
one must resort to approximations or to computer sim- 
ulations. Unfortunately, the latter only give access to 
free energy derivatives and free energy differences. Sev- 
eral advanced techniques have been developed that allow 
to relate free energies of different state points to each 
other, and a large body of literature has been written 
on this topio^'^'^'^'^'ii^. Nevertheless, comparing the free 
energies of arbitrary systems remains a challenge, and 
alternative approaches that allow to determine the abso- 
lute free energy for each individual system are clearly of 
interest. 

On principle, absolute free energies can be obtained by 
connecting the system of interest with a reference system 
of known free energy. In this letter, we propose a general 
strategy for the construction of analytically solvable ref- 
erence systems, that can be connected with a wide class 
of structures via thermodynamic integration. 

Thermodynamic integration^iii^ is a widely applied 
method to determine free energy differences. The ba- 
sic idea is the following: Consider a system of N parti- 
cles with a Hamiltonian H{r^ ,e), which explicitly 
depends on some parameter e. In order to obtain an ex- 
pression for the free energy of the system, one uses the 
relation dF/de = {dH{e) / de) , where (...) denotes the 
thermodynamic average. Here and in the following, we 



set /csT — 1. In general, {dH{e) / de) , is directly acces- 
sible in a simulation. Thus the expression above can be 
used to evaluate the free energy difference between two 
systems at different e: One samples {dH{e) / de) for a 
range of e and integrates 

AF^F{e,)-Fieo)= I" de'l^^^y . (1) 

If the free energy is known for one £o (reference system) , 
the method can be used to calculate absolute free ener- 
gies for a whole range of e. However, it is crucial that 
the evolution of {dH{e) / de) on the integration path is 
reversible, i.e., no phase transition of first order may be 
crossed. This limits the choice of suitable integration 
paths and reference systems. For gases the ideal gas is a 
useful reference system, for crystals the "Einstein crys- 
tal" (a crystal where the particles are bound to sites of 
a fixed lattice by harmonic springsiiii^) . To the best 
of our knowledge, no general reference system has been 
introduced so far that can be used for arbitrary dense 
disordered systems. 

Our central idea to remedy this situation is very sim- 
ple. We propose to take a configuration that is represen- 
tative for the structure of interest (obtained e.g. within a 
typical simulation of an equilibrated system) and to con- 
struct a reference system by first 'pinning' this configura- 
tion with suitable external fields, and then switching off 
the internal interactions. In the remainder of this letter, 
we will show how this idea can be exploited to evaluate 
absolute free energies in practice. 

For the purpose of illustration, we begin by considering 
the Ising model i/o = ^JJ2{ij) ^i^j, where {ij) denotes 
neighbouring i and j and Si = ±1. To evaluate the free 
energy Fq at a given temperature, we simulate the system 
until it is equilibrated, and then pick one typical config- 
uration {sf^} as 'representative' reference configuration. 
The reference system is then defined by the Hamiltonian 

H,,f{e)^-eJ2s^s^ , (2) 

i 

and its free energy can be computed easily, Frei{e) = 
— A^ln (2 cosh (e)). To establish the connection with 
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the original system, we procede in two steps: First we 
define an intermediate model H'{e) — Hq + H,-f,{{£), 
which reduces to Hq at e = 0. The free energy differ- 
ence AFi(e) — Fq — F'{e) between the original system 
and the intermediate system can be calculated for arbi- 
trary e by thermodynamic integration, using {dH' /de) = 
— {^iSisf). We choose e large enough that the spins 
in the system H'{e) hardly fluctuate about the reference 
value. The second step is to connect the intermediate sys- 
tem to the reference system. The free energy difference 
between the two systems at the same value of e, AF2 (e) = 
— Fref (e), is evaluated by carrying out a simulation 
with additional Monte Carlo (MC) moves that switch 
on and off the interaction J according to a Metropolis 
criterion. We obtain AF2{e) = — ln(Pon/^off), where 
^on,off is the fraction of configurations with interactions 
switched on (rsp. off) in the simulation. Combining ev- 
erything, we can finally calculate the absolute free energy 
of the target system Hq, Fq ^ Frcf{e) + AFi{e) + AF2{e). 

Now we transfer this idea to off-lattice particle mod- 
els. For clarity, we only discuss monatomic liquids and 
solids in the NVT ensemble in the following. Our method 
can easily be generalized to molecular systems, and, as 
we shall demonstrate below, to constant pressure simula- 
tions. Furthermore, we disregard the kinetic contribution 
to the free energy, which can be evaluated trivially 

Let configurations be characterized by a set of coordi- 
nates {r,;} and the configurational energy be given by a 
Hamiltonian Hq = U{{ri}). To calculate the free energy 
of a given, arbitrary equilibrium structure, we choose a 
'representative' configuration {rf"}, obtained, e.g., from 
a simulation of an equilibrated system, and construct a 
reference system by imposing local potentials 
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that pin the particles' positions to the reference posi- 
tions r^. Here <& defines attractive potential- wells cen- 
tered at each position rf^, with ^{x) < for a; < 1 
and $ = elsewhere. Note that particle i can only be 
trapped by well i and not by the other wells. To make the 
particles indistinguishable, as they should be, we allow 
them to swap identities (i.e., labels i,j) at regular in- 
tervals during the simulations. We will show below that 
such identity swaps are also necessary to equilibrate the 
system efhciently. 

The (configurational) reference free energy is given by 
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where Vq is the volume of the sphere of radius rcutoff and 
(7$(e) := d dx x'^~^ (e"'^*^'^) — l) for a d-dimensional 
problem. In our simulations, we mostly used a linear well 
potential, {x) = x — 1. In this case one has (7$(£) = 

d/e'^ ^e^ — X]fe=o ^'^ 1^^^ ■ As before, we also define an in- 
termediate model H'{e) — Ho + H,-c{{e), and evaluate the 
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Figure 1: Sketch of moves in our Monte Carlo algorithm, (a) 
Simple particle displacements. (Could be replaced e.g., by 
short Molecular Dynamics runs.) (b) Smart particle swaps, 
(c) Smart particle relocations. See text for explanation 



free energy difference between the reference system and 
the intermediate system at high e with a MC simulation 
where the interaction Hq is switched on and off (if neces- 
sary, in several steps) . The free energy difference between 
the target system and the intermediate system is com- 
puted by sampling dF'/de = - ''/^l/'^cutoff)) 
and performing a thermodynamic integration. The re- 
maining challenge is to devise an algorithm for sampling 
the intermediate model efficiently for arbitrary e. 

Before describing such an algorithm, we briefly com- 
ment on the relation between our method and the 
Einstein crystal method to determine free energies of 
solids^. In the Einstein crystal method, the particles 
are not swapped, and the reference system is a regu- 
lar lattice of harmonic wells with infinite cutoff radius. 
This works well as long as the particles in the target sys- 
tem stay close to their respective well positions. In a 
liquid, where their mean-square displacement diverges, 
{J2i ^ilfi — r^l/rcutos)) diverges as well for small e and 
can no longer be sampled. Therefore the introduction of 
a finite cutoff is crucial. We note that our method can 
also be used to evaluate free energies of crystals. 

Setting a finite range for the reference potential, how- 
ever, introduces a different problem: The particles need 
to find their respective wells of attraction. We therefore 
introduce two MC moves that help particles i explore 
their well i (Fig.[T|). One move (Fig. [l]b) swaps particles 
in a smart way. It works as follows: 

• Pick a random particle i and find the set of particles 
{n-i} that are within the attraction range of well i. 

• If particle i ^ {ui}: pick a particle j from {rii}, and 
swap i and j with the probability min{l, j^e~'^^ }. 

• Otherwise: pick a particle j from all particles 



if j ^ {rii}: swap with probability min{l, — 
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- if J G {rii}: swap with probability min{l, e^^^ }. 
Here AH' is the difference of the energies (according to 
the intermediate model) of the old and new configura- 
tion. This algorithm promotes particle swaps that bring 
particles close to their respective well and nevertheless 
satisfies detailed balance. 

The other move (Fig. [T] c) relocates particles i with a 
bias towards the neighborhood of their well i: 

• Pick a random particle i (with position r^). 

• Choose a new position from a given (biased) distri- 
bution P,{r'i) = exp{-W{\r'^ - rf |)). 



3 




Monte Carlo Steps 



Figure 2: Illustration of the effect of the moves of Fig. [T]on 
the equilibration of a system of 80 hard disks (diameter D) 
at a density p — 0.8/ D^, after switching on linear well po- 
tentials with strength e = 50 (rcutoff ~ 2D). Swap moves 
and relocation moves (one per bead) were attempted one per 
100 MC sweeps. Left: Evolution of {<E>) in simulations that 
include different moves as indicated. Right: Corresponding fi- 
nal configurations. Circles indicate particle positions, crosses 
give well positions. Particles and their respective wells are 
connected by straight lines. 



• Relocate the particle from r to r'^ vifith probability 

min{l,P(ri)/P(rOe-^^'}. 
Obvious choices for W{r) which we have tested are 
W{r) = £$(r/rcutoff), or W{r) = const, for r < rcutoff- 
At high £, the relocation move helps to overcome trapped 
situations where most particles are bound to a well, and 
a few cannot escape from a local cage. To illustrate the 
effect of the different moves, Fig. [5] shows the evolution 
of the observable averaged over all particles i, in a 
two dimensional system of hard disks, after e had been 
raised from zero to a high value. In a MC simulation that 
includes only random particle displacements, the system 
is far from equilibration after one million MC sweeps (a). 
The smart swap moves speed up equilibration consid- 
erably, but the system gets trapped in a configuration 
where one particle cannot enter its well (b). This prob- 
lem is solved by including smart relocation moves (c). 

Table I: Results for the free energy of hard spheres. F/N^^ 
is the value according to the Carnahan-Starling equation of 
stated, a) linear potential $, liquid reference state, b) linear 
^, hep reference state, c) harmonic liquid reference state. 



N/V F/N_ (F/iV)cs 

0.25 0.62(0) 0.625 

0.5"' 1.54(1) 1.544 

0.5*"' 1.54(0) 1.544 

0.5^=' 1.54(9) 1.544 

0.75 3.00(9) 3.005 



We will now demonstrate the power of our approach at 



a few examples. They are not meant to be self-sufRcient 
scientific studies - the simulations were carried out on 
simple workstations with poorly optimized test programs 
in relatively short time (less than a week). We are aware 
that a careful analysis of the finite-size effects should be 
done in all cases. Here we only intend to illustrate some 
properties of the algorithm. 

We have studied hard spheres in two (2d) and three 
dimensions (3d). For the remainder of this letter we use 
the particle diameter D as unit of length. Table U shows 
results for the free energy of a 3d liquid of hard spheres. 
The simulations were performed on a system of iV = 256 
particles, using 50 values of e and 6 ■ 10^ MC sweeps for 
each value at N/V = 0.25 and N/V = 0.5, and 200 val- 
ues of s times 1 Mio. MC sweeps at N/V = 0.75. The 
results agree with the values obtained by integration of 
the Carnahan-Starling equation of stated within the er- 
ror bars. For N/V = 0.5 we compared the cases a) linear 
potential $ and liquid reference state, b) linear $ and 
crystalline reference state, and c) harmonic $ and liquid 
reference state. Within the error-bars these variations 
produce the same result. However, for more accurate 
calculations the linear potential seems to be most useful, 
because the particles get trapped most efficiently. In case 
b) we did not see a hysteresis on increasing/decreasing 
e. Apparently, the trapping of particles in a crystalline 
array of wells is not associated with a phase transition at 
this density. This will presumably be different closer to 
liquid/solid coexistence. Nevertheless we can conclude 
that our method is quite robust and may work even if 
the reference configuration is not 'ideal', i.e., not repre- 
sentative of the target structure. 

Next we show an example for the application of the 
method to dense disordered systems, where the dynam- 
ics is driven by cooperative processes. We studied hard 
disks in 2d up to densities where the equilibrium phase 
is a solid, and enforced a vacancy defect by taking one 
particle out of an otherwise ordered configuration. These 
simulations were carried out at constant pressure P in a 
rectangular simulation box of varying area, but fixed side 
ratio 1 : \/3/4, to accomodate a triangular lattice. The 
defect is then stable, but highly mobile. 

We compare three different structures (Fig.[3|): An or- 
dered solid (a), an ordered solid with a vacancy (b), and a 
metastable disordered jammed phase (c), which was ob- 
tained by compressing the system from the fluid phase. 
Free energy calculations were carried out at P = 10 for 
these three cases, and additionally at P = 6, in the fluid 
regime. To calculate the free energy in the enthalpic en- 
semble, we use a reference system that is deflned in terms 
of scalable coordinates {i.e., the positions of the well cen- 
ters are rescaled along with the particle coordinates if the 
volume of the system changes), and pin the volume of the 
system by an additional term e{V — l^of)^ in the refer- 
ence Hamiltonian. The thermodynamic integration was 
carried out for e = ... 40 using 66 values of e, and 10 
Mio. MC sweeps for each e. The resulting free enthalpies 
G can be related to the chemical potential fi by virtue of 
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Density [1/D ] 



Figure 3: Characterization of the dense two-dimensional sys- 
tems discussed in the text. Left: Pressure vs. density as 
obtained from constant pressure simulations at e = 0. (a) 
A'^ = 100 particles, expanded from an ordered solid phase, (b) 
A'' = 99 particles, expanded from an ordered solid phase with 
one vacancy, (c) N = 100 particles (diameter D), compressed 
from the fluid phase. The solid line shows the theoretical 
estimate P = p/(l - ^pf with p = {N + 1)/{V). Right: 
A configuration with one vacancy at the beginning (crosses) 
and the end (circles) of a MC run. The thin and thick arrows 
mark the position of the defect at the beginning and the end. 



the thermodynamic relation G — /iiV. 

To set the frame, we show in Fig. [3] a) the pressure- 
density curves for the cases a),b), and c). In the fluid 
regime (p < 0.8), they can be fitted nicely with the theo- 
retical estimatei'^ P = p/(l-f p)^ with p ^ {N + 1)/{V). 
Fig. Ob) illustrates the mobility of the defect at pressure 
P = 10. It should be noted that in the solid regime, 
the center of mass motion of the complete system is not 
sampled well, because individual particles stay close to 
their lattice sites. A similar problem is encountered in 
the Einstein crystal method and has lead to the devel- 
opment of the 'Einstein molecule—, where the reference 
crystal is defined in terms of relative coordinates. This 
idea can easily be transfered to our approach. Here, we 
ignore the center-of-mass correction, because it is smaller 
than our statistical error. 

At P = 6, the free energy calculation yields the free 



enthalpy per particle fj, = 8.997(2), which is reason- 
ably close to the theoretical estimate fi — 9.047 ob- 
tained by integrating the theoretical equation of state. 
At P = 10, we found /isoUd — 13.617(2) in the solid 
state, and pjam — 13.675(2) in the jammed state, which 
establishes that the solid is indeed the stable phase. For 
the system with one defect, we obtained the total en- 
thalpy Gdcfoct = 1361.7(2). This result can be used 
to estimate the core free energy of the vacancy fic — 
Gdcfcct — ^J'so\idN + lii{N) = 7.1(3), which corresponds to 
a relative vacancy frequency of roughly 10""^. (For com- 
parison, the frequency of vacancies at liquid/solid coex- 
istence in 3d is roughly 10~4ii.) Probably Hc is largely 
overestimated due to finite size effects, hence the value 
given above is at best an upper bound. More detailed 
studies shall be carried out in the future. Here, the ex- 
ample mainly serves to illustrate the use of our approach 
in situations where free energies are difficult to access 
with other methods. 

In summary, we have introduced a general method to 
compute absolute free energies for a wide range of struc- 
tures. We have illustrated the method for monatomic 
simple systems, but we believe that it can be applied 
equally well to molecular fiuids and mixtures. We an- 
ticipate that our method will be useful to calculate free 
energies of systems that are not directly connected with 
the ideal gas, such as liquid crystal phases, membranes, 
or proteins in solution. Also defects in solids seem to 
be a promising field of application. From a fundamental 
point of view, it should be interesting to study how well 
the method can be applied to glassy systems, which have 
not just one, but a whole set of 'representative' config- 
urations, one for each local minimum in the rugged free 
energy landscape. 
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